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Abstract. We study quasi-static deformation of dense granular packings. In 
the reference configuration, a granular material is under confining stress (pre- 
stress). Then the packing is deformed by imposing external boundary condi- 
tions, which model engineering experiments such as shear and compression. 
The deformation is assumed to preserve the local structure of neighbors for 
each particle, which is a realistic assumption for highly compacted packings 
driven by small boundary displacements. We propose a two-dimensional net- 
work model of such deformations. The model takes into account elastic inter- 
particle interactions and incorporates geometric impenetrability constraints. 
The effects of friction are neglected. In our model, a granular packing is rep- 
resented by a spring-lattice network, whereby the particle centers correspond 
to vertices of the network, and interparticle contacts correspond to the edges. 
We work with general network geometries: periodicity is not assumed. For 
the springs, we use a quadratic elastic energy function. Combined with the 
linearized impenetrability constraints, this function provides a regularization 
of the hard-sphere potential for small displacements. 

When the network deforms, each spring either preserves its length (this 
corresponds to a solid-like contact), or expands (this represents a broken con- 
tact). Our goal is to study distribution of solid-like contacts in the energy- 
minimizing configuration. We prove that under certain geometric conditions 
on the network, there are at least two non-stretched springs attached to each 
node, which means that every particle has at least two solid-like contacts. The 
result implies that a particle cannot loose contact with all of its neighbors. 
This eliminates micro-avalanches as a mechanism for structural weakening in 
small shear deformation. 



Key words: granular materials, constrained optimization, geometric rigidity, dis- 
crete variational inequalities. 

1. Introduction 

Materials that are composed of collections of separate, macroscopic solid grains 
belong to the general classification of granular materials. Examples of such ma- 
terials are common, including sand, gravel, medicinal pills, coins, and breakfast 
cereal. Granular media are important to numerous industries ranging from min- 
ing to pharmaceuticals. In geophysics, granular materials are a central problem in 
understanding the physics of earthquakes and tectonic faulting. Earthquake fault 
zones produce granular wear material continuously as a function of shear and grind- 
ing between the fault surfaces. The wear material, known as fault gouge, varies in 
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thickness from 10's of cm to 1000 m and plays a critical role in determining the 
fault zone frictional strength, the stability of fault slip, and the size of the rupture 
nucleation dimension. 

Granular media display a variety of complex static and dynamic properties that 
distinguish them from conventional solids and liquids. The complexity of granular 
media lies primarily in the collective properties of a macroscopic number of grains 
and how they interact with each other. The conditions under which a granular 
medium is stable or flows and the nature of this flow depend critically on the 
distributions of grain size and shape as well as the interactions between the grains. 
The practical importance of granular media combined with the richness of their 
physical properties has led to a great deal of interest from both theoretical and 
experimental points of view |H 1111 IT2*] . 

An important class of granular materials consists of nearly rigid particles that 
possess the following property: if a moderate force is applied, the particles start 
to move, and only after a substantial increase of the force the particles deform 
significantly. In other words, for loads which are not very high, the deformations 
inside a particle are small compared with the displacement of the particle center 
of mass. Consequently, the particle shapes change very little, so that each particle 
can be associated with a region of space that is inaccessible to any other particle. 
This gives rise to constraints on the admissible positions of particles. These impen- 
etrability constraints are also known as geometric, kinematic, and excluded volume 
constraints. 

A physical phenomenon related to appearance of constraints is jamming. A par- 
ticle is jammed when its motion is completely obstructed by the neighbors, so the 
whole cluster of neighboring particles can only move together as a rigid body. The 
corresponding mathematical notion of rigidity ( 23] ) can be applied to various phys- 
ical (sphere packings, frameworks (trusses)), as well as mathematical objects. In 
particular, an important mathematical object associated with any particle packing 
is a contact graph defined as follows: vertices of this graph are particle centers of 
mass, while edges represent interparticle contacts. 

The simplest physical model that exhibits jamming is a classical hard sphere 
packing. The particles in this model are represented by rigid spheres, and the only 
interparticle forces are reactions of constraints. Rigidity of hard sphere packings 
is studied in 0. This problem can be formulated as a problem of detecting rigid- 
ity of the cable framework associated with the contact graph of the packing. The 
framework is obtained by replacing edges of the contact graph with the cables, 
and vertices with flexible hinges. The lengths of the cables can increase but not 
decrease, which models the impenetrability constraints. Recently, a linear program- 
ming algorithm for detecting rigidity in hard sphere packings (equivalently, cable 
frameworks) was proposed in [5]. 

In this work, we also use bar frameworks. A bar framework is obtained from 
a graph by replacing the edges with rigid bars, and vertices with hinges. A bar 
framework and the associated graph are called rigid if the only possible vertex 
motions correspond to rigid body motions of the whole framework. We note that 
both bar and cable frameworks can be associated to the same graph. To generate 
the bar framework, the edges of a graph are replaced with rigid bars that can only 
translate and rotate. In the case of the cable framework, one replaces edges with 
cables that can either move as rigid bodies or stretch. Thus every motion of a bar 
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framework is also permitted by the cable framework, but the converse is not true 
in general. Therefore, it is possible that a bar framework associated to a graph 
is rigid, while the cable framework corresponding to the same graph is not. Both 
cable frameworks and bar frameworks are special case of the so-called tensegrity 
frameworks studied in In a tensegrity framework, properties of edges can vary, 
e.g. some edges may be bars, others may be cables or struts (that can shrink, but 
not stretch). 

It appears that the currently available mathematical results |2J [31 |S] on statics 
of discrete particle systems with geometric constraints deal only with hard particle 
packings. To the best of authors' knowledge, there are no results on frictional pack- 
ings, and even elastic frictionless packings have yet not been studied. The present 
work differs from El IS] m several respects. First, all these studies deal with 
rigid particles. We consider a somewhat more realistic situation of geometrically 
constrained particles with elastic interactions defined by a quadratic potential en- 
ergy. Second, while 0|21E] focus on jamming, we are interested in generic contact 
patterns of the the energy minimizing configurations. The packings that we study 
are not jammed. Their contact graphs are such that the associated bar framework 
is rigid, but the packing can still deform when external boundary conditions are 
applied. 

The third difference is in the type of the boundary conditions. The conditions in 
[5] are periodic or hard wall conditions. The periodic conditions are commonly used 
to minimize influence of the boundaries in the problem. However, presence of walls 
is a major factor that determines bulk behavior of granular materials. Therefore it 
seems better to use boundary conditions corresponding to engineering and physical 
experiments, where the walls are rigid and may be moving. 

A frequently observed property of granular materials is concentration of the bulk 
deformation in thin layers called shear bands. Within a band, the contact forces are 
weak, and the relative displacements can be on the order of particle size or larger. 
For quasi-static flows driven by small shear rates, the corresponding patterns are 
called micro-bands f|15j). In that paper, shear band structures were studied by 
means of numerical simulations. The simulations in |15| show that the typical size 
and number of bands in quasi-static shear depend on the imposed shear rate. For 
small shear rates, the bands have length and width comparable with the particle 
size. The distribution of these micro-bands within the material is rather uniform. 
As the shear rate increases, the band structure exhibits coarsening: the number 
of bands becomes smaller, and the length of each band increases. For sufficiently 
large shear rates, a single macroscopic shear band appears. 

Here, we are interested in the case of small external boundary conditions. In 
that case, a micro-band can be formed by weakening of a single contact, or a 
small group of neighboring contacts. All weak contacts form a subnetwork of the 
whole contact network. Such networks of weak contacts (corresponding to the 
micro-band patterns in |15j) were studied numerically in 20 . The goal of this 
paper is to describe some generic geometric features of micro-band, or weak contact 
networks in dense packings of nearly rigid particles. In small deformation, pattern 
formation may be caused by the local jamming (which mathematically amounts 
to impenetrability constraints), and friction. We are concerned with the role of 
constraints, while friction in neglected. The notion of high density at this point is 
rather intuitive. It could mean, for instance, that each particle is in contact with at 
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least three other particles, and the packing is jammed in the reference configuration, 
subject to zero boundary conditions. Below we make the notion of high density 
more precise (see the second paragraph on p. 4), using the relationship between 
the contact graph and the Delaunay graph (see e.g. |7j), generated by the set of 
particle centers. 

In two dimensions, particles are represented by disks Di of radii with centers 
x l , i = 1,2,..., TV. The initial reference configuration is deformed by applying 
prescribed small displacements to the boundary particles. Assuming that the de- 
formations inside of the individual particles are small, and neglecting rotational 
degrees of freedom, one can characterize the deformations of Di by the displace- 
ments u l of their centers. The elastic interaction forces are modeled as in classical 
mechanics of point particles: the force exerted by Dj on Di is applied at x 4 , its 
direction is along the line joining x l and x J , and its magnitude depends linearly on 
u J — u l . 

We further assume that the granular material is pre-stressed (or, equivalently, 
the material is under confining stress). This means that in the reference state, 
the particles are squashed into each other as a result of applied external pressure. 
Further compression is supposed to be impossible (requires infinite energy), which 
introduces impenetrability constraints into the problem. To model impenetrability, 
one can, for instance, require that 

(1.1) |(x i + u l )-(x^+u^')|> a . 1 + aj , 

for each pair of particles. Since the the packing is dense, and the displacements 
are expected to be small, it makes sense to require that a particle cannot escape 
a cage formed by its neighbors. Therefore, the contacts that exist in the reference 
configuration may be broken, but no new contacts are created after applying exter- 
nal boundary conditions to the reference configuration. An important consequence 
of this assumption is as follows. If the displacements satisfy (0.1) for each pair of 
particles in contact, then for any pair of particles (0.1) is automatically satisfied. 
Indeed, if two particles are not in contact in the reference configuration, they cannot 
come into contact in the deformed configuration, and the distance between them 
must be larger than the sum of their radii. In the sequel, we use this assumption 
in the course of proving the main result of this paper ( Theorem 15. Q . 

Next, we introduce a network model which describes a granular material under 
the above assumptions. The vertices of the network are the particle centers, and the 
edges represent particle contacts. The collection of vertices x ! , i = 1, 2, . . . , N, and 
edges forms the contact network (graph) T. We suppose that T is a triangulation 
of a connected, convex polygonal domain f2. This assumption is realistic, since, for 
example, a periodic 2D packing of disks is triangular. Another natural triangulation 
generated by x% i = 1,2, ... ,N is the Delaunay graph G. In principle, T and G 
may be different, since some edges in G may not correspond to contacts. In the 
present case, we suppose that T and G coincide, which corresponds to "maximally 
dense" packings. 

For small displacements u 1 , the quadratic constraints (|1.1|) can be approximated 
by their linearizations near u 1 = u 7 =0, which leads to the linearized impenetra- 
bility constraints 



(1.2) 



ill' u:q >0, t = l,2,...,N 
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for each pair of vertices i,j connected by an edge of T. In (|1.2|l . q u = (x- 7 — 
x 4 )/ jx- 7 — x 4 | are unit vectors that point from x* to x- 7 along the line of centers. 
Note that if the position of Di is fixed (u J = 0) , then u J satisfying 1)1.2(1 must lie 
in the half-plane u • q 2 - 7 > 0, so that Dj would be moving away from Di. 

For certain boundary conditions the deformed packing can become more loose 
than the reference packing. On the macroscale, this can be observed as swelling of 
the specimen caused by the increase in the volume of the void space between the 
particles. Such swelling is typical in shear deformation, where the overall volume 
increase, known as dilatation, is observed in experiments. To increase the void 
volume, some of the contacts present in the reference configuration must be broken 
in the deformed configuration. Therefore, among all the contacts (satisfying (jl^ L 
we further distinguish two types of contacts: broken and solid-like (see Fig. 1). We 
call a contact broken if 

(1.3) {u j - u*) • q ij > 0, 
and solid-like if 

(1.4) (u j - u*) • q lJ = 0. 

The solid-like contacts correspond to two possible types of pair motions. The 
first type is a rigid motion of a pair, in which case the contact is called stuck. 

-e 



Figure 1. Left: a broken contact. Center: a solid-like sheared 
contact. Right: a solid like stuck contact (infinitesimal rigid rota- 
tion of a pair). The arrows indicate displacements. 

The second type is a local shear motion. In the local coordinates of one particle, 
it is either the motion of the second particle in the direction perpendicular to the 
line of centers, or an infinitesimal rotation (rolling). The corresponding contacts are 
called sheared. In our idealized model, friction in neglected, and any tangential force 
would lead to immediate separation of particles, because for disk-shaped particles, 
the contact surface is a point. We, however, still call these contacts solid-like, 
because in reality, these contacts are subject to friction forces, the contact surface 
has a positive area, and the particles in a sheared contact will not separate until 
the tangential force reaches the static friction threshold. In the case of rolling, the 
particles stay in contact and the pair is capable of bearing a compressive load. 

Physically, vertices of the network can be realized as unit point masses and edges 
can be realized as elastic springs. Elastic force of the spring is determined by 
the pair potential H(tij), where = (u- 7 — u 1 ) • q*- 7 . The potential is an important 
ingredient of our model, and therefore we discuss it in detail. To motivate the 
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choice of H, we first recall the classical hard sphere potential Hhs, which in our 
notation is defined by 



(1.5) H hs {U 



oo if tij < 0, 
if t i:j > 0. 



Hhs models the following two options: (i) moving non-deformable (hard spheres) 
particles toward each other requires infinite energy (a vertical line at Uj = 0), (ii) 
moving particles apart requires no energy. Note that (|1.5f) already incorporates the 
constraints 1|1.2[) by requiring infinite potential energy to violate the constraints. 




Figure 2 . a) Hard sphere potential Hhs \ b) The elastic potential 
H(tij) combines a vertical wall at tij — and a quadratic function 
with the vertex at d. 

Elastic interaction between Di and Dj, together with constraints (|1.2|) can be 
modeled by the following potential 

Here d characterizes the cut-off distance of the potential, and C determines the 
magnitude of the pre-stress potential (the value of the potential when tij = 0. 
The potential l|1.6|) is shown on Fig. 2, together with the hard-sphere potential. 
The formula describes two options: (i) moving particles toward each other 

requires infinite energy; (ii) movement of particles apart from each other is caused 
by finite, linear elastic force f y = —dH/dii 1 ; This force is repulsive for small 
distances (tij < d), since dH/dt 1 ^ < 0. The magnitude of the force / <J '(d) = 
dH ^' d ~> = C^d~ 3 \tij — d\. The magnitude of the force of pre-stress (or confining 

stress) is given by lim t . ,^ + / y = ^Cd~ 2 , which tends to zero as d — > oo and C 
is fixed. Therfore, the effect of pre-stress is smaller for larger d. Further, H(Uj,d) 
regularizes Hh s in the following sense: if d > da > 0, then lim^— ►«> H (tij, d) = 
uniformly on (0, do]. In the paper, we do not pass to this limit. Instead we choose d 
sufficiently large and fix it, so that H is close to Hhs- Also, for technical simplicity, 
we set C — 1 in i|1.6|) . which corresponds to an appropriate rescaling. 

In reality, once the distance between Di and Dj is greater than the sum of their 
radii a, + dj, the pair interaction force is zero. In our model, we still have a small 
repulsive force for all ai+a,j < Uj < d. So, the particles in our model would continue 
accelerating away from each other even after separating. This favors separation 
of particles, and could lead to increase in the number of broken contacts. Since 
our goal is estimating the number of solid-like contacts from below, this increase 
is acceptable. In fact, our estimate holds for all d sufficiently large. We also 
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mention that elastic contact force predicted by the classical Hertz theory is a non- 
linear function of tjj . Our model is chosen for simplicity, and can be viewed as an 
approximation of Hertz theory, valid for sufficiently small displacements. 

In general, the cut-off parameter d may be different for different pairs of particles 
in contact. Therefore, we define a pair interaction energy 

(1.7) KU^dij) = -d-^Uj-dij) 2 , 

and consider 

„, , , f 00 if tin < 0, 

(L8) HiuM = [ h(tiM [{ - , ( ; 

where 

(1.9) dy ; = <%, % € [1/2,1]. 

The formula (f 1 . H|l is more general than (|1.6(1 . The choice of dy in 1)1. 9JI ensures 
that for all pairs (i,j), the points where h(tij,dij) — and located in the interval 
[\/2d,d\. The number 1/2 is of no particular significance. Any number in the 
interval (0, 1) would work just as well. We only need all to have the same sign 
and comparable magnitudes controlled by d. Finally, the total elastic interaction 
energy Q of the network is obtained by summing up h(Uj, dij) over all pairs (i, j) 
corresponding to the edges of the network. 

The equilibrium state of a granular material corresponds to a minimum of Q 
subject to the constraints l|1.2|) and the appropriate boundary conditions. Since the 
functional Q is quadratic, and the constraints and boundary conditions are linear, 
this is a quadratic programming problem, studied in optimization theory (e.g., |10|). 
In the language of optimization theory, the solid-like contacts (|1.4|) correspond to 
the so-called active constraints, while the constraints corresponding to the broken 
contacts (|1.3|l are called inactive. The question addressed in this paper concerns 
the total number and spatial distribution of each type of constraints in the energy- 
minimizing configuration of the network. It appears that no general results of this 
type are currently available in optimization theory. The present study makes use of 
the geometric features of the contact graph, in particular its rigidity properties (see 
above), to investigate the energy minimizer. We also note the connection between 
our constrained variational problem and continuum variational inequalities [5], 
|14|. |17|. Our problem can be viewed as a discrete variational inequality. 

The main result of the paper is Theorem 15.11 in Section El There we consider 
a packing whose contact graph in the reference configuration is a triangulation of 
a convex polygonal domain. The packing is deformed by imposing displacement 
boundary conditions at the packing boundary. The boundary conditions model the 
motion of rigid walls in engineering experiments. We prove that for generic contact 
graphs (precisely defined in Definition 14. 21) and generic pre-stresses (corresponding 
to the choice of Sij in 1)1. 7J) . the constrained energy minimizer for sufficiently large 
d provides a packing with at least two solid-like contacts per each particle. There 
also some (non-generic) choices of 5{j for which theorem does not provide a definite 
conclusion. 

The network of solid-like contacts is the load-bearing structure. The network of 
broken contacts can be associated with the micro-bands |lf>) , |16| that appear during 
small shear deformations. The result implies that no particle can lose contact with 
all of its neighbors, which eliminates "micro-avalanches" . Put another way, loss of 
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structural integrity in dense packings is evolutionary rather than catastrophic, so 
that shearing with a small displacement will first lead to dilatation, during which 
the packing becomes more loose everywhere, and only then local avalanches may 
occur. 

Another useful consequence of theorem 15.11 is as follows. It provides a lower 
bound on the order parameter, recently introduced in P3 I22| as one of the main 
ingredients of the new phenomenological theory of dense granular flows proposed by 
Aronson and Tsimring. The order parameter p characterizes the phase transition 
from solid to fluidized state. To define p at an arbitrary point x of Q, one begins 
by fixing a mesoscopic averaging volume V of characteristic size h (e.g., a disk of 
radius h centered at x). Then, all solid- like contacts within V should be counted. 
Next, the obtained number n s of solid- like contacts is divided by the number n of 
all contacts within V to obtain p. So, p is a mesoscopic average, which in general 
depends on h. In many systems, such as periodic elastic composites, the results 
of mesoscopic averaging is practically independent of h for h larger than a certain 
characteristic length. For disordered granular materials this is not necessarily true. 
Therefore, a rigorous mathematical theory may require study of a family of order 
parameters parametrized by h. In Section 6, we define such a family of order 
parameters using the notion of a fc-neighborhood of a vertex of T as a discrete 
analogue of V. Specifying an integer k in our definition corresponds to choosing h 
in the continuous case. 

The paper is organized as follows. In Sect. 2 we formulate the main constrained 
minimization problem. The problem contains two types of constraints: impenetra- 
bility constraints (|1.2|l and the boundary constraints (see H2.6|) . Ij2.7|l ). correspond- 
ing to the external boundary conditions. Elimination of these boundary constraints 
leads to a reduced minimization problem. In Sect. 3 we recall some facts concerning 
first-order rigidity of graphs. In Sect. 4 we show existence of a unique minimizcr 
of the reduced problem. Optimality conditions for the reduced problem are stated 
and analyzed in Sect. 5, where we also state and prove the main theorem 15. II In 
Sect. 6 we introduce a definition of the order parameter in the spirit of |22| and 
give a lower bound on the order parameter that follows from the main theorem. 
Finally, conclusions are provided in Sect. 7. 



2. Formulation of the problem 

2.1. Elastic interactions with impenetrability constraints. In 2D, consider 
a packing of spheres Di of radii a% with centers x l , i — 1, 2, . . . , N. (All vectors in 
this paper are column vectors and we use superscript 'T' to indicate transposition). 
The packing fills a bounded region. After an infinitesimal motion, the position of 
the center of Di is y\ We write y l = x l + u 1 where u l are displacements. The 
vertices x 1 , x J are connected by an edge if and only if Di, Dj are in contact. In this 
case, we call x' and x- 7 neighbors. We denote by Mi the set of j e {1,2,..., N} 
such that x J is a neighbor of x\ Orientation of contacts (equivalently, edges) is 
prescribed by the unit vectors 

(2-1) q« = 

\xP — x l | 

The vertices x* and edges (i, j) define the contact graph T. Let E denote the number 
of edges of T. The edge set £ of V is given by {(£, j) : j G M i} i = 1, 2, . . . , N}. To 
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each edge we can associate a pair potential energy h(Uj,dij) defined in l|1.7jl . 
Summing up all these energies we obtain the total elastic interaction energy of the 
network. It is a quadratic form 
(2.2) 

N N 

q(u\ u 2 , . . . u n ) = y Y h (tij,dij) = -d- 3 y E (( uJ - ui ) • ^ ~ d ^) 2 > 

i=l j£Af; i=l jGAfi 

on the displacements u l , i = 1, 2, . . . , N. In (|2.2|) . d, rfy are parameters specified by 

(Ol) . nn . 

Our objective is to determine the displacements u 1 , i - - 1, 2, . . . , N so that the 
the energy functional Q is minimized subject to two types of constraints. The 
first type of constraints consists of linearized impenetrability constraints. These 
are obtained by formally linearizing the condition that the distance between two 
spheres in contact cannot decrease. Consider two spheres Di, Dj in contact. In the 
reference configuration, 

(2.3) |x' -x J '| = a, t +aj. 
Assuming that Di, Dj cannot overlap, we have 

(2.4) \r -yj\ > ai + aj . 

These are the impenetrability constraints. We linearize ({2.4(1 by writing 

ly'-y^'l 2 = |x J -x J +u l -u j \ 2 

= |x l - x j \ 2 + 2(x i - x j ) ■ (u* - u J ) + |u l - u j \ 2 
= (a t + aj) 2 + 2(x s: - x J ) • (u J - u J ) + |u l - u J | 2 . 

Now for for "small" |u l — u J | we can neglect quadratic term |u l — u J | 2 , and (|2.4|l 
yields 2(x l — x- 7 ) • (u l — u J ) > 0, which in turn is equivalent to (u J — u*) • q y > 
where q u is as defined in 12.1|l . Therefore, the first set of constraints we impose on 
the displacements u 1 , u 2 , . . . , is 

(2.5) (u j - u J ) • > 0, j G Mi, i = 1, 2, . . . ,N. 

The second type of constraints corresponds to the boundary conditions. Particles 
located at the packing boundary have prescribed displacements. In the sequel we 
refer to these particles as boundary particles. The corresponding vertices of T are 
called boundary vertices. Other particles are referred to as interior, or sometimes, 
free, and the corresponding vertices of T as interior vertices. 

All boundary particles are divided into several groups, numbered 1,2,... , M. 
Let I m denote the set of indices of the particles in group m for m = 1, 2, . . . , M. 
Each sphere in a certain group is in contact with at least one other sphere from 
the same group. Each group moves as a single rigid body. We assume that the 
prescribed boundary displacements are of the form 

(2.6) u l = R m (x 4 ), i G 7 m , m = 1, 2, . . . , M, 
where 

(2.7) R m (x l ) = c m + a m K (x l - x*< m ), i G I m , m = 1, 2, . . . , M, 



10 K. A. ARIYAWANSA, LEONID BERLYAND, AND ALEXANDER PANCHENKO 



and c m , x* ,m are given vectors, a m is a given scalar, and K is the matrix denoting 
clockwise rotation by 7r/2. The functions R m are called infinitesimal rigid displace- 
ments, parametrized by a scalar a m , and vectors c m and x*' m . We refer the reader 
to Sect. |21 for more details on rigid displacements. 
Our description above leads to the 



Main problem: 

(2.8) minimize ^(u 1 , u 2 , . . . , u N ) 

(2.9) subject to linearized impenetrability constraints l|2.5(l 

(2.10) and boundary conditions H2.6[l . 

2.2. Feasible region. Let us define the configuration space U. Points of this space 
are denoted by U = ((u 1 ) 1 ", (u 2 ) T , . . . , (u Ar ) T ) T . 

Remark. To avoid this heavy notation, we simply write 

U = (u 1 ,u 2 ,...,u JV ), 

when no confusion can occur. 

Dimension of U is 2N. Feasible region T is the subset of U in which all the con- 
straints H2.5fl and l|2.ti[) are satisfied. The points satisfying (|2.5|l form a polyhedral 
(not necessarily bounded) region. The boundary of this region consists of parts of 
the hyperplanes (subspaces of dimension 2N — 1) defined by 

(2.11) (u* - u*) • q y =0, jeAfi, * = 1,2,...,JV. 

Because of the close relation to rigidity, we refer to l|2.11|l as R- equations. Equations 
(12. 6|) define M planes S m , m — 1, . . . , M. Dimensions of S m depend on the number 
of the boundary particles in the m-th group. 

For each point of U G jF, some of the constraints l|2.5[) are satisfied as equations. 
These constraints are called active. The corresponding edges of the contact graph F 
are called active as well. The rest of (|2.5|l are satisfied as strict inequalities. These 
are inactive constraints (respectively, edges). 

2.3. Elimination of constraints corresponding to boundary conditions. 

The quadratic form Q in (|2.2() can be written in a convenient form in terms of a 
certain matrix R r . To define R r , we index the edges oi T by I, I — 1,2, .... E. Let 
S £ be the edge of T corresponding to I for I = 1, 2, . . . , E. Let R r be the 
E x 2N matrix whose /-th row is defined by 

+(q«<)i if m = 2(ji - 1) + 1 

+(q i «') 2 if m = 2(7i-l) + 2 

-{cC dl )i ifm = 2(ij-l) + l 

-(q iljl ) 2 if m = 2(ii - 1) +2 

otherwise 



(2.12) R\ n 



for 1 = 1,2,..., E. 

Remarks. 1. R r is the (first-order) rigidity matrix, a well known object in geo- 
metric rigidity theory (see e.g. (3"ll2"3*p. 

2. Consider vertices x 2! ,x j! and the edge I connecting them. The corresponding 
row r l of R r has 2N entries. We can view v l as a string of N pairs of numbers, the 
first pair corresponding to x 1 , the second to x 2 and so on. For simplicity, we shall 
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call a pair of entries corresponding to a particular vertex x 1 a place corresponding 
to x\ 

Then we can interpret equation (|2.12|) as follows. A row r l has zeros at all 
places, except two. The non-zero entries are — q* iJi , written as a two-dimensional 
row vector at the place corresponding to x l! ; and q 1 ' Ji , written as a two-dimensional 
row at the place corresponding to x Jl . 

3. A row of R r corresponds to an edge of V. Therefore it is natural to call a 
row active (respectively, inactive) if a corresponding edge is active (respectively, 
inactive) . 

Now define the vector d € M. E by 

(2.13) d = —(di 1 j l: di 2 j 2: di E j E ), 

where di u j t are chosen according to ljl.9|l . With these notations the quadratic form 
Q in (|2.2|1 can be written as 

(2.14) Q(U) =d- 3 ^(R r V + d) • (iTU + d). 

We now eliminate the boundary conditions i|2.6[l from the main problem (|2. 8112. 91 
I2.10fl . Let Nb — J2m=i card(I m ). Then the equations l|2.6|l simply state that the 
2Nb components of U corresponding to the Nb boundary vertices have prescribed 
displacements. Without loss of generality assume that the last 2Nb components of 
U correspond to the boundary vertices. Let us partition U as 

(2.15) U=[z|w], 

where z = (u 1 , u 2 , . . . , u N ~ Nh ) corresponds to displacement vectors of interior ver- 
tices, and w = (u N ~ Nb+1 , i+2 , , , . ; U N ^ corresponds to the displacements of 
the boundary vertices. The equality constraint (|2.6|l is now simply 

(2.16) w = g 

where g € M. 2Nb is the vector of displacements prescribed by the right-hand-sides 
of (|2.6|l . The matrix R r can be partitioned similarly to 12.15f) : 

(2.17) R r = [ R | R h ], 

where dimensions of R and R b are E x 2(N — Nb) and E x 2Nb, respectively. Denote 

(2.18) a = i? b g. 

Using l|2.15|l - (|2.18|) in (|2.14|) and in l|2.5|l we can reduce the main problem (|2. 81 12.91 
l2~TUjl to 

Reduced problem: 

(2.19) minimize F (z) = ^d~ 3 (Rz + a + d) • (Rz + a + d) 

(2.20) subject to Rz + a > 0. 

The minimization in (|2.19|l is taken over all z G M.( N ~ Nb *> . 
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3. First-order rigidity 
A rigid motion is a composition of a translation and rotation: 
(3.1) y(x) = c + x* + 0(x-x*), 

where O is an orthogonal (rotation) matrix, c is a translation vector, x* is a center 
of rotation. If O is close to identity I (infinitesimally small rotation), then 

O ~I + A, 

where A is a skew matrix (a^ = —dji). 

Suppose that in a two-dimensional rigid motion, the rotation angle a is close to 
zero. Then 

0=( C ° Sa Sina W I ")+a( ° l)=I + aK, 
V — sma cos a / V 1 / V — 10/ 



where 



K 



1 

-1 



is a clockwise rotation by 7r/2. In that case, ll.'i. l[l becomes 



(3.2) y(x) = c + x + aK(x 



c\ + x* — a(x - 
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c 2 + x\ + a(x — x*)i 
Let u = y(x) — x denote the displacement. We can write H3.2JI as 

(3.3) u(x) = c + aK(x-x*). 

Definition 3.1. We call (|3.3|l an infinitesimal rigid displacements in 2D. 

Next, let G be a graph. Consider all motions of vertices of G that preserve the 
lengths of the edges. If the only such motions are the rigid body motions of the 
whole graph, then the graph is called rigid. A graph T is first-order rigid 23 if all 
solutions of the i?-system (|2.11ll are infinitesimally rigid displacements. In addition, 
r is independent if the rows of the rigidity matrix R r are linearly independent. 
Graphs that are both first-order rigid and independent are called isostatic f|23|). 
Intuitively, an isostatic graph is minimally rigid, that is removing any edge results 
in loss of rigidity. Another notion of rigidity is generic rigidity, see |23) . According 
to Thm. 49.1.7 from [22], generic rigidity for a neighborhood in a configuration 
space is equivalent to the first-order rigidity for some specific configuration in that 
neighborhood. 

With the definition of R r and U in Sect. Othe system (|2~TT|) can be written 

as 

(3.4) R r U = 0. 

Note the connection between i?-system and constraints, as well as the functional 
of the main problem $M .9lll0|) . 

The following definition (see |23| for a d-dimensional definition) is useful for 
verifying rigidity of graphs. 

Definition 3.2. For a graph T, the Henneberg 2-construction in 2D is a sequence 
of graphs G\,G%,...,G n such that: 

(i) Gfc+i is obtained from Gk by either vertex addition (attaching a new vertex by 
2 edges); or edge splitting (replacing and edge from Gk with a new vertex joined 
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to its ends and to 1 other vertex); 

(ii) Gk is a complete graph on k vertices, and G n = Y . 

The following result is stated in ([22], thm. 49.1.13): 

Theorem 3.3. If a graph G cl 2 is obtained by a Henneberg 2 -construction, then 
G is generically isostatic. 

Remark 1. In veiw of the definitions above, Theorem 13.31 implies that a graph ob- 
tained by Henneberg 2-construction is first-order rigid and independent (minimally 
rigid). 

In the present case, the rows of the rigidity matrix R r are not linearly indepen- 
dent, but the row rank is maximal. This means that we typically have more edges 
than needed to ensure rigidity of T. In this situation, the following theorem (thm. 
49.1.14 from [23]) is useful. 

Theorem 3.4. If two graphs G\ and Gi are generically rigid planar graphs sharing 
at least 2 vertices, then the graph G obtained by combining all vertices and edges of 
G\ , G2 is generically rigid. 

Remark 2. Because of the relation between generic rigidity and first-order rigid- 
ity, Theorem 13.41 implies that combining first-order rigid graphs G\, G2 as in this 
Theorem yields a first-order rigid graph G. We shall use Theorem 13.41 to obtain 
first-order rigidity of triangulations. Indeed, one triangle G\ is first-order rigid. 
Adding another triangle G2 so that G\ and Gi share an edge, yields a first-order 
rigid graph. Then we can proceed sequentially. Given a first-order rigid triangu- 
lation Gfc, we construct Gk+i by combining Gk with a triangle. This new triangle 
either shares two vertices with Gk, or all three vertices. In the first case, Gk+i 
would have one more vertex and two more edges than Gk- In the second case, 
Gk+i has the same number of vertices as Gk, and one more edge. By theorem 13. 41 
any planar triangulation obtained by this sequential procedure is first-order rigid. 

4. Existence and uniqueness of minimizers of the reduced problem 

Let Q be a bounded connected domain in R 2 with a polygonal boundary. First 
we show that, under certain assumptions on geometry of T, the matrix R has full 
column rank. 

We shall say that Y is a triangulation if edges of Y partition £1 into a disjoint 
union of triangles. 

Let X be a set of interior vertices of Y, containing at least two elements. Consider 
a graph Gx C Y defined as follows. Vertices of Gx are all elements of X. Edges of 
Gx are those edges of Y that join two vertices from X . We also assume that X is 
chosen so that Gx is a connected graph. 

Definition 4.1. The contact graph Y is cell- connected if for each Gx C Y as above, 
there exist two vertices x-^x 2 in Gx, and two interior vertices x^x 2 in Y\Gx, such 
that the quadrilateral with vertices x 1 , x 2 , x 1 , x 2 is a union of two adjacent triangles 

of r. 

An example illustrating the definition if shown in Fig. 3. 
Definition 4.2. We call Y a regular triangulation if 

(i) r can be obtained by sequential addition of triangles in the way described in 
the Remark 2; 
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Figure 3. a) -A cell-connected graph. Solid dots indicate vertices 
of Gx ; b) - The cell connection with three edges and four vertices 
is shown separately; c)- A triangulation that is not cell connected. 
Here the vertices of the small triangle inside are connected to other 
vertices by pairs of collinear edges. 



(ii) every interior vertex, connected to a boundary vertex, is also connected to at 
least one other boundary vertex, and the corresponding edges are non-collinear; 

(iii) T is cell-connected. 

Remark. Note that i) in Definition 14.21 implies that Y is first-order rigid. The 
property ii) states that every edge connecting a boundary vertex with an interior 
vertex must be a part of the boundary of a triangle, containing two boundary 
vertices and one interior vertex. 

Informally, Definition ^, ll f or (iii) in Definition ^. 21 ) can be interpreted as a strong 
connectivity property. According to Definition 14.11 a generic connected subgraph 
Gx is connected to the "rest of Y" not just by an edge, but by a "more robust" cell 
structure that consists of three edges, with one edge bracing the other two, (see Fig 
3.). Note also that either x 1 or x 2 are connected to the vertices of Y\Gx by a pair 
of non-collinear edges. This observation is important in the proof of Proposition 
14.31 below. It is not difficult to see that a periodic triangular planar graph satisfies 
iii). However, there are triangulations with a mean coordination number four that 
do not satisfy iii). An example of such a graph is shown in Fig. 3 c). 

Proposition 4.3. Suppose that Y is a regular triangulation. Then rank R = 
2(N-N b ). 

Proof. Consider a subgraph Y max C Y constructed inductively as follows. Begin 
with Ti that consists of all boundary vertices. On the next step, add an interior 
vertices connected to Ti by two or more non-collinear edges. Also, add exactly two 
non-collinear edges that connect this vertex to Ti. Call the resulting graph T2. 
Generally, given T^, k > 2, define Yk+i — Yk U Sk, where Sk consists of an interior 
vertex x fc , not contained in Yk but connected to Yk by at least two non-collinear 
edges, together with a pair of non-collinear edges connecting x fc to Yk- Since the 
graph r has a finite number of vertices, the process terminates after a finite number 
of steps. The resulting graph is Y max . The construction is illustrated in Fig. We 
claim that Y max contains all vertices of Y. To obtain a contradiction, suppose that 
there are vertices not included into Y max . Denote the set of these vertices by X, 
and denote by Gx the graph formed by vertices in X and all edges of Y that connect 
these vertices. Let G c x be any connected component of Gx- If G c x is a single point 
x s , then, since T is a triangulation, there must be at least three edges incident at 
x 9 , and at least two of these edges must be non-collinear. Thus x g must be in Y max , 
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FIGURE 4. The subgraphs Ti (left) and T 2 (right). The boundary 
vertices are shown in gray. The interior vertices of Ti and T2 are 
shown in black. The edges of Ti, T2 are shown by solid lines. 
The other vertices of F are represented by the unfilled circles. The 
edges of T not included into Fi,r2 are represented by the dotted 
lines. 





Figure 5. Two different subgraphs T max constructed inductively. 
Both subgraphs are constructed starting with T2 in Fig. 4. 



which gives a contradiction. Next suppose that G c x contains two or more vertices. 
By Definition 14. II (see also the Remark following that Definition), there must be a 
pair of vertices of x 1 , x 2 in G c x connected to two vertices in T\Gx by three edges, 
and at least one of x*,x 2 must be connected to r\Gjs: by two non-collinear edges. 
Denote this vertex by x s . It must be included into T max which gives a contradiction 
and proves the claim. 

Next, we claim that the number of edges in T max is 2(N — Nb). Indeed, each 
interior vertex in T2 has exactly two non-collinear edges incident at it. Then, on 
each of the next steps, we add an interior vertex together with two non-collinear 
edges incident at it. Since the number of free vertices in T max is N ~ Nb, the claim 
is proved. 

Finally, we claim that the rows of R corresponding to the edges of T max are 
linearly independent. Let the matrix of these rows be denoted by R m ax- This is 
a square 2{N — Nb) matrix. We claim that an appropriate row-reduction reduces 
Rmax to a matrix R' max that has block-diagonal form: for each vertex x l there are 
exactly two rows r i,:l ,r 1 ' 2 of R' max , and two linearly independent unit vectors q yi 
and q ya , such that r 1,1 (r 1,2 ) contains q Ul (q lj2 ) at a place corresponding to x 1 
while all other entries in these rows are zero. 
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To see this, consider first a "basic unit" of T2: an interior vertex x 1 and two non- 
collinear edges incident at it. Let the corresponding unit vectors be q 11 , q* 2 . Recall 
that these edges connect x l to two boundary vertices. Consequently, the rows r 1 , r 2 
corresponding to the above pair of edges have zeros at all places, except two places 
corresponding to x\ The non-zero entries of r 1 (r 2 ) are two components of q l1 
(q l2 ). Since q 11 , q* 2 are linearly independent, so are r x ,r 2 . Furthermore, linear 
combinations of r 1 ,! -2 can be used to eliminate non-zero entries in other rows. By 
adding an appropriate linear combination of r 1 , r 2 to a row with some unit vector 
q y at a place corresponding to x 1 , we can obtain zeros at this place. Hence, by using 
r 1 ,! -2 as pivots in Gaussian elimination we can obtain rows whose only non-zero 
entries are at a place corresponding to a vertex that was added to T2 on the next 
step of the iteration. These rows, in turn, can be used as pivots. Continuing with 
row reduction, we can eventually reduce all rows of R m ax to this form (this follows 
from the fact that T max contains all vertices of T). The proposition is proved. 
Remark. It is interesting to compare R and the rigidity matrix R r . It is well- 
known that for a first-order rigid graph, the null space of R r is non-trivial and 
consists of infinitesimal rigid displacements. The proposition above shows that 
the null space of R is trivial. The main difference in structure between these two 
matrices is that R contains special rows, that might be called broken. These rows 
correspond to edges connecting an interior vertex to a boundary one. A typical 
row of R r has four non-zero entries, while each broken row has only two. These 
entries occur at a place corresponding to an interior vertex. If an interior vertex is 
connected to two boundary vertices, then the regular triangulation property of T 
ensures that the corresponding broken rows are non-collinear, and can be used in 
the sequential Gaussian elimination, as done in the proof. 

Proposition 4.4. Consider the problem \2. 1912.20)) . Suppose that T is a regu- 
lar triangulation, the feasible set of \2.19\2.2(]\) is non-empty, and that the uncon- 
strained minimizer of F(z) is not feasible. Then the problem $2.iyi2.2(l\) admits a 
unique minimizer that is a point on the boundary of its feasible set. 

Proof. The problem (|lli)li>.lWjl has a feasible point z. Then the problem 

has a unique minimizer z* as we now demonstrate. Let the set C(z) be defined by 

£(z) = {z £ M 2 ^-^) : F(z) < F(z)}. 

By Proposition 14.31 the matrix R has full column rank. Therefore, the set £(z) is 
an ellipsoid, which is a closed, bounded, convex set. The set TL of z 6 R 2 ^-^) 
satisfying the constraint (|2.2(J|) is a closed half space. Therefore, £(z) n TL is a 
nonempty, closed, bounded convex set. Indeed, the reduced problem 12.1912.20(1 is 
equivalent to the problem 

(4.1) minimize F(z) 

(4.2) subject to z e £(z) n TL. 

Now by the continuity of F and the compactness of £(z) n TL we see that the 
problem I|4.1I4.2|I and hence the reduced problem I|2.19I2.20|I has a minimizer z*. 
Since R has full column rank, R T R is positive definite. The positive definiteness 
of the matrix R T R implies that F is strictly convex on £(z) n TL., from which we 
conclude that z* must be unique. 

Now if i?z*+a > 0, then z* must be the unconstrained minimizer of F. This con- 
tradicts the assumption that the unconstrained minimizer is not feasible. Therefore, 
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some components of Rz* + a must be zero, and thus z* must be on the boundary 
of the feasible region. 

5. OPTIMALITY CONDITIONS FOR THE REDUCED PROBLEM 

The reduced problem (|2.19I2.20|) has a convex objective function and linear con- 
straints. For such problems it is possible to state optimality conditions that are 
both necessary and sufficient ^B]. To be specific define the Lagrangian L for the 
problem I|2.19I2.20[) by 

(5.1) L(z, A) = -cT^ite + a + d) ■ (Rz + a + d) - eT 3 A • {Rz + a). 

Then z* solves problem (|2.19I2.2U|I if and only if there exists A* such that z = z* 
and A = A* satisfy the Karush, Kuhn, Tucker (KKT) conditions 

V z £(z, A) = 
Rz + a > 
A • (Rz + a) = 
A > 

or equivalently 



(5.2) R (Rz + a + d- X) = 

(5.3) i?z + a>0 

(5.4) A • (Rz + a) = 

(5.5) A>0. 



See [El Chapter 12]. 

Before stating and proving the main result (Theorem 15.11) , we list all the as- 
sumptions, including both new and previously used. 

Al. Consider the problem of minimizing (|2.2(l subject only to the boundary condi- 
tions H2.6II2.7|) . but not the constraints (|1.2|) . We assume that the minimizer of that 
problem is not feasible, that is, this minimizer does not satisfy the impenetrability 
constraints (|1.2f) . Further we assume that the feasible region is not empty, which 
means that there is at lest one point z satisfying all the inequality constraints l|5.3|l . 
A2. The network T is a regular triangulation (as defined in Definition 14.2(1 . 
A3. The boundary conditions are prescribed so that 

(5.6) |(x l + u ! ) - (x* + u j )\ < a,i + a 7 - + min a k , 

J fe=l,2,...iV 

for each pair x 4 , x- 7 of boundary vertices in contact. 

Let us provide some comments on the nature of assumptions A1-A3. Assump- 
tion Al means that minimizing the energy of the spring network subject only to 
boundary conditions leads to a configuration in which at least one spring is com- 
pressed (and thus violates the impenetrability constraints). 

Assumption A2 concerns the contact geometry. The edges of the network split 
the domain of the problem (a polygon) into elementary cells (triangles). Near the 
boundary, the cells must be compatible with the geometry of the boundary in the 
following sense. If an interior vertex is connected to a boundary vertex, then it is 
also connected with another boundary vertex, located next to the first boundary 
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vertex. Hence, every triangular cell adjacent to the exterior boundary must contain 
one free vertex and two boundary vertices. 

Assumption A3 means that the boundary conditions l|2.6l 12.7( 1 are chosen to 
prevent particles from escaping through the gaps made by displacing the boundary 
particles. Clearly, if two boundary particles belong to the same group, then no gap 
can appear between them, and |(x* + u*) — (x 1 + u J )| = a t + ctj. Formation of gaps 
would be possible between two boundary particles from different groups which are 
in contact in the reference configuration. If the parameters of rigid body motions in 
the boundary conditions 1(2.61 12.71) are prescribed arbitrarily, then the two particles 
may move away from each other, and open a gap large enough for a third particle 
to slip through. Assumption A3 prohibits formation of such gaps. 

Theorem 5.1. Suppose that assumptions A1-A3 hold. Then there exist d* > 
and a vector S = (#i, 82, ■ ■ ■ , Se) £ K B : Si G (— 1, —1/2), I = 1, 2, . . . ,E, such that 
for each d = dS with d > d* , the unique minimizer of l2.lfJl2.2Cfy has the following 
property. Each interior vertex x 1 of T has at least two active edges incident at it. 
The corresponding unit vectors q 1 '- 71 , q Ij2 must be linearly independent. 

Proof. 

Step 1. We claim that A3 implies that there is cq > 0, which depends on the 
boundary conditions, but is independent of the choice of in 12.2J1 . such that 
each feasible displacement u 1 ,i = 1,2, . . . , N, satisfies 

(5.7) K\<co, 

k = 1,2. Indeed, first we observe that if u 1 ,^ satisfy the linearized constraint 
1)1.2(1 . then they also satisfy the distance constraint (|l.l(l (the converse is not true 
in general) . Then any feasible collection of displacements also satisfies the distance 
constraints 1 (1.1(1 for each pair of neighboring vertices. Now we recall the assumption 
made in the introduction to conclude that l|l.lfl must hold for all pairs of vertices. 
Fix I G {1,2,..., N}, corresponding to an interior vertex, and consider a smaller 
packing V of particles, containing only Di and all boundary particles. In the ref- 
erence configuration, Di is completely surrounded by boundary particles. Then 
the boundary conditions are prescribed according to A3, the boundary particles 
still completely confine Di , so that x' must displace to x' + w' that lies inside a 
certain bounded domain fl 1 that depends only on boundary conditions. Since xj. 
are bounded, this implies that the claim is true for all displacements w' which are 
feasible for the smaller packing V . Clearly the set of all such displacements is larger 
than the set of all u' feasible under all constraints (|l.l(l , and the latter set is larger 
than the set of all u' feasible under the linearized constraints (|1.2(l . This proves 
the claim. 
Step 2. Let 

(5.8) v* = ]T q ij , i=l,2,...,2(N-N b ). 

First, we prove the theorem under the additional assumption 

(5.9) For each i = 1, 2, . . . , 2(N - N b ), v l ^ sq lj , 

where j £ A/j, s e E. 

We note that l|5.9|l implies that 

(5.10) |v*| > v > 
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with vo independent of i. Indeed, s can be zero, so validity of (|5.9(l means in 
particular that all v l are non-zero. Since there is finitely many v l , 1|5.1()JI holds. 

Consider solutions of the KKT system (|5.2I5.3I5.4I5.5|) . From (|5.4[) . (|5.5|) it 
follows that Xj = if the j-th constraint is inactive. Let Oj — (Rz + a)j. If 
the j-th. constraint is active then Oj — 0, while Xj is arbitrary. Suppose that 
a feasible point z* is given. Then Oj are given. To solve (|5.2|) we need to find 
A. Denote by r k ,k = 1,2,..., E the rows of R (the columns of R T ), and sup- 
pose that the rows r , r 2 ...,r correspond to the active constraints, and that 
the rows r s+1 , r s+2 , . . . , r E correspond to the inactive constraints. Choose d = 
(— d, —d, . . . , — d). Then (|5.2|) can be written as 

S E E 

(5.11) -£VAi+ J2 r'ft + d^r' = 0. 

1=1 l=S+l 1=1 

Pick a vertex x l of T and consider the restriction of each r' in (|5.11(1 to the two 
components corresponding to x l . Then we have 

(5.12) - W + W + dvl = °' 

where the first sum is taken over active edges incident at x% while the second sum 
is over the inactive edges incident at x\ 

Next, we determine the minimal number of active edges needed for H5.12(l to hold. 
We can look at (|5.12J) as a local problem in which u 1 may vary, while u J , j G M are 
fixed. Denote by J-i C K 2 the feasible region of this local problem. By A3, Ti is a 
polygon, each side of which corresponds to one or more constraints being active. 

In the generic case, one constraint per side is active. In the non-generic case, 
two or more active constraints correspond to the same side. Since our goal is 
estimating the number of active constraints from below, it is sufficient to consider 
only the generic case, corresponding to the "worst case scenario". In the generic 
case there are only three possibilities. 

Case 1. u 1 is inside Ti. All edges incident at x l are inactive. 

Case 2. u l belongs to only one of the sides of dTi. One edge is active. 

Case 3. u l is a vertex of J-i. Two edges are active. 

Consider case 1. Then Ij5.12ll cannot hold for d sufficiently large. Indeed, |v l | > 
t;o > by assumption, while | X^sa/' ^ij^l 1S bounded from above independent of 
d in view of l|5.7(l . 

Consider case 2. Let us number the active edge by (i, 1). Then (|5.12(l can be 
written as 

(5.13) -Aaq' 1 + - u ') ' I* 3 ') 1 W + = °' 

Enlarging d, if necessary, we see that (|5.13|l can hold only if 

(5.14) v^sq^ 1 , 

where s < 0. Since Ij5.14|l is not allowed by H5.8|l . (|5.12|l cannot hold for sufficiently 
large d. 
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Consider case 3. Number the two active edges by (i,2). The equation 

lETTl is 

(5.15) -A a q l1 - A i2 q 12 + ('"' u : q")q'' • r/v' (I. 

isA/< ,i>2 

For this to hold for large d, v l must be a non-positive linear combination of q l1 , q l2 . 
These two vectors are linearly independent, otherwise their intersection would not 
be a vertex of Ti. So, Case 3 is possible, provided v l lies in the negative cone of 
two active edges. 

Step 3. Now we remove the assumption (|5.9|l . For each i — 1,2, ... , (N — Nb), and 
each 8 G R E , define a two-dimensional vector v l to be the restriction of R 8 = 
12f=i <^ r ' to a place i. The theorem will be proved is we show that there is a choice 
of 8 such that 8i G [1/2, 1], and v l has property (|5.9(l . Indeed, if such 8 is found, 
we could choose d = dS, where d > is sufficiently large, and repeat the arguments 
made in the first step, using v l instead of v l . 

To show existence of 8, consider the cube Ce = {y G M b : yi G (1/2, 1),/ = 
1,2, ... ,E}. Pick any point y* G Ce- Since Ce is open, there is a Euclidean 
open ball B(y*,p) C Ce, with the radius p > 0. Consider the image of B(y*,p) 
under the mapping R . Since R has full rank, R is surjective, and is therefore 
an open mapping. Thus, R (B(y* , p)) contains a Euclidean open ball B(R y* , p*) 

of a positive radius p* depending only on R and p, but not on y*. If R y* has 
property (|5.9|) . we choose 8 = y* and we are done. Otherwise, note that for each 

i = l,2,..., (N — Nb), the ball B(R y* , p*) contains a non-empty two-dimensional 

Euclidean open ball Bi centered at the restriction of R y* to the place i. Since 
for each i the set {v G R 2 : v = sq lJ , s G M.,j G A/" 1 } is a union of a finite 
number of lines, it cannot contain a two-dimensional ball. Therefore, for each 
i = 1,2,..., (N — Nb) there must be a vector v l G Bi having property (|5.9|l . Now 
we can define YliLi ^/ r ' y i a its restrictions v\ Next, by construction, we can find 

T e 

a vector S G B(y*, p) c Ce such that R 8 = J2i=i ■ 

The theorem is proved. 
Remark. The choice of 8 in the proof is based on the following criterion. Consider 

the vector R T 8 G W. 2( - N ~ Nb \ The proof works if 

(5.16) R T 8 ^ (vi, v 2 , . . . , Vjv-jvJ, 

where at least one of the two-dimensional vectors v^, i = 1,2,..., (N — Nb) is of the 
form V; = sq JJ , for some real s and j G Ni- In other words, if is inadmissible, 
then it must lie on the line through the origin with direction vector q y . For each 
fixed i, the inadmissible set Vi = {v^ G M 2 : = sq u } has Hausdorff dimension one 
(a finite union of lines in the plane), while the admissible set is the two-dimensional 
complement of Vi in R 2 . Therefore, the set of inadmissible vectors in the right hand 
side of 1)5.16(1 has dimension 2 (TV — Nb) — 1 while the set of admissible vectors R 8 

is of dimension 2(N — Nb). Thus the admissible R 8 are generic, and the theorem 
holds for a generic (in the sense of Definition 14. 2(1 packing under a generic pre-stress 
(the latter is determined by a generic choice of 8). 



a network model for granular statics with impenetrability constraints 

6. Order parameter 

Recently, a phenomenological theory of slow dense granular flows was proposed 
in [Tll22|. A key quantity in that theory is the order parameter, defined as the ratio 
of the number of solid-like contacts to the number of all contacts within a given 
control volume. In 22 , a contact is considered solid-like if two particles are jammed 
together for longer than a characteristic collision time. The relevant characteristic 
time is r = a/v a , where a is particle radius and v a is the speed of sound in a solid 
material of the particles. Our model corresponds to the instantaneous material 
response, when r is much smaller than other relaxation times in the system, such 
as the ratio of the sample size to a typical particle velocity. 

An obvious type of pair motion leading to a solid-like contact is a rigid dis- 
placement (a pair of particles infinitesimally moves as a rigid body). We shall 
call this type of contacts stuck. If a contact between Dj and Dj is stuck, then 
(u 1 — u J ) • q y = 0, which is easy to check using the definition of rigid displace- 
ments. This means that the impenetrability constraint for the corresponding edge 
of of the network is satisfied as an equation (the edge is active). However, not 
every active edge corresponds to a stuck contact. Another type of a local motion 
that produces (u l — u J ) • q y = is an infinitesimal shear motion when u 1 — u J is 
orthogonal to q u . The corresponding contact is called sheared. Note also that infin- 
itesimal shear is the same as infinitesimal rotation, so this type of motion includes 
infinitesimal rolling as well as shear sliding. 

We consider both sheared and stuck contact as solid-like, because stuck contacts 
are stable, while sheared contacts in an actual granular material will be subject to 
friction. Friction can be viewed as partially stabilizing, at least when the shear- 
ing force is below the static friction threshold. Such non-sliding frictional contacts 
are considered as solid-like in the simulations performed in [22} • In addition, some 
heuristic arguments and numerical simulations presented in [3 EI j suggest that fric- 
tion enhances elastic behavior of sufficiently large samples. Therefore, it makes 
sense to think of the network of solid-like contacts as the main load-bearing struc- 
ture and call this network strong. In contrast, a broken contact satisfying 
corresponds to a local weakening in the material because in this case two particles 
separate completely. We can think of the network of all broken contacts as weak. 
Moreover, division of contacts into broken and solid-like corresponds to the division 
of constraints into active and inactive, as done in optimization theory. Therefore, 
this division is natural mathematically, and also makes sense from the physics point 
of view. 

In addition, the definition in [22] does not sufficiently clarify the nature of av- 
eraging. The notion of an order parameter in static problems should not use time 
averaging. The result of spatial averaging depends on the size of the sample that is 
being averaged. Thus, if the order parameter is obtained by, say, spatial averaging, 
then it must depend on both location and size of the "control volume" . In the dis- 
crete situation, the size of the averaging sample can be measured by the minimal 
number of edges connecting a pair of vertices within the sample. 

This suggests a definition of the size-dependent order parameter. To state this 
definition we first define the averaging sample. 

Definition 6.1. A vertex x- 5 is in the k-th neighborhood of x 1 if T contains a path 
connecting x 1 and x J with no more than k edges. 



22 K. A. ARIYAWANSA, LEONID BERLYAND, AND ALEXANDER PANCHENKO 



Now, to each /c-neighborhood we can associate a value of an order parameter. 

Definition 6.2. For each x* and each non-negative integer k < N, the size depen- 
dent order parameter p(x l , k) is defined by 



(6.1) p(x\fc) 



E* 



where the numerator is the number of active edges in neighborhood of x 1 , and 
denominator is the number of all edges in that neighborhood. 

Theorem 15 . II implies the lower bound 

N 

(6-2) p(x\A0>- 

on the order parameter associated with the maximal, iV-th neighborhood of each 
interior vertex x\ Indeed, counting active edges (two per vertex) gives 2N edges, 
each counted at most twice. In particular, Ij6.2|l means that the order parameter 
p(x J , N) is bounded from below by the reciprocal of the mean coordination number 
of the network. 



7. Conclusions 

We have studied a network model of quasi-static deformation of dense pre- 
stressed granular materials. In our model, the packing was represented by a network 
of linear elastic springs. Each spring corresponds to a contact between two par- 
ticles. Geometric impenetrability constraints within the packing were modeled by 
the linearized impenetrability constraints on the displacements of the vertices of 
the network. The constraints have the form of linear inequalities, that can be sat- 
isfied either as an equality (an active constraint) , or as a strict inequality (inactive 
constraint). Constraints are in one-to one correspondence with the interpaticle 
contacts. An active constraint corresponds to a relatively stable solid-like contact. 
Inactive constraints represent the relatively weak broken contacts. The question 
addressed in the paper is to estimate the total number and distribution of the 
solid-like contacts in the energy-minimizing configuration. We showed that each 
interior vertex of the network has at least two solid-like contacts corresponding 
to it. This result qualitatively reproduces the micro-band structure obtained in 
[151 116| by numerical simulations. We also discussed the connection between our 
result and a lower bound on the order parameter [Tll22|. In the paper, we proposed 
a definition of the order parameter that is similar to the one introduced in |22j . 
but differs from it in the interpretation of the so-called solid-like contacts. On the 
one hand, our definition appears to be in accord with a physical picture of granular 
statics, recently proposed in [HIE]. On the other hand, it is a naturally related to 
optimization theory. 
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